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To understand the observational properties of cosmological models, in particular, the temper- 
ature of the cosmic microwave background radiation, it is necessary to study their null geodesies. 
Dynamical systems theory, in conjunction with the orthonormal frame approach, has proved to be 
an invaluable tool for analyzing spatially homogeneous cosmologies. It is thus natural to use such 
techniques to study the geodesies of these models. We therefore augment the Einstein field equa- 
tions with the geodesic equations, all written in dimensionless form, obtaining an extended system 
of first-order ordinary differential equations that simultaneously describes the evolution of the grav- 
itational field and the behavior of the associated geodesies. It is shown that the extended system 
is a powerful tool for investigating the effect of spacetime anisotropies on the temperature of the 
cosmic microwave background radiation, and that it can also be used for studying geodesic chaos. 



1. INTRODUCTION 

The dynamical systems approach to the field equations of general relativity has been an invaluable tool for gaining 
qualitative information about the solution space of the anisotropic but spatially homogeneous (SH) Bianchi cosmologies 
(see Wainwright & EUisQ j^] and references therein). In this approach one uses the orthonormal frame formalism of Ellis 
& MacCallum |^ to write the field equations as an autonomous system of first-order differential equations, the evolution 
equations for the gravitational field. One can then apply techniques from the theory of dynamical systems to obtain 
qualitative information about the evolution of Bianchi cosmologies. The essential step is to introduce dimensionless 
variables for the gravitational field by normalizing with the rate-of-expansion scalar, or equivalently, the Hubble 
scalar. A consequence of this choice of variables is that the equilibrium points of the evolution equations correspond 
to self-similar Bianchi models, leading to the insight that this special subclass of models plays a fundamental role 
in determining the structure of the general solution space. An added bonus is that the evolution equations are well 
suited for doing numerical simulations of Bianchi cosmologies. 

In order to understand the observational properties of the Bianchi models, however, it is necessary to study the 
behavior of their null geodesies. In this paper we augment the evolution equations of the gravitational field with 
the geodesic equations using the components of the tangent vector field as the basic variables, thereby creating an 
extended system of equations. This yields a system of coupled first-order ordinary differential equations that describes 
the evolution of the gravitational field and the behavior of the associated geodesies. It turns out that normalizing 
the geodesic variables with the energy leads to bounded variables for null and timelike geodesies, which is of great 
advantage. 

It is widely believed that a highly isotropic cosmic microwave background (CMB) temperature implies that the uni- 
verse as a whole must be highly isotropic about our position, and thus accurately described by a Friedmann-Lemaitre 
(FL) model. Bianchi cosmologies provide an arena for testing this belief. Since the 1960s, various investigations of 
the CMB temperature in SH universes have used the observed anisotropy in the temperature to place restrictions on 
the overall anisotropy of the expansion of the universe, as described by the dimensionless scalar]^ E = cr/{y/3H) (see. 



^From now on we will refer to this reference as WE. 

■^Here — ^OabO'"''^ is the norm of the shear tensor and H is the Hubble variable. 
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for example, Collins & Hawking Q). Some of these investigations have also determined the temperature patterns on 
the celestial sphere in universes of different Bianchi types (see, for example. Barrow et al. [Q). The studies that have 
been performed to date, however, suffer from a number of limitations: 

i) They are restricted to those Bianchi group types that are admitted by the FL models. Indeed, the most detailed 
analyses, for example, Bajtlik et al. [p|, have considered only the simplest Bianchi types, namely, I and V. 

ii) The results are derived using linear perturbations of the FL models. Such a simple approach cannot be justified 
in all situations (see Collins & Hawking page 316 and Doroshkevich et al. [gj, page 558). 

iii) The analyses provide no bounds on the intrinsic anisotropy in the gravitational field, as described, for example, 
by a dimensionless scalar W formed from the Weyl curvature tensor (see Wainwright et al. 0, page 2580, for 
the definition of W). 

The extended system of equations is a powerful tool for investigating the anisotropy of the CMB temperature free of 
the above limitations. In particular, the method can be applied even if the model in question is not close to an FL 
model. 

The outline of the paper is as follows: In section || we show how to extend the orthonormal frame formalism to 
include the geodesic equations in SH Bianchi cosmologies. As examples we consider diagonal class A models and type 
V and type VI;i models of class B. In section || the structure of the extended system of equations is discussed. Section 
^ contains examples of the dynamics of geodesies in some self-similar cosmological models. As a simple non-self-similar 
example we consider the locally rotationally symmetric (LRS) Bianchi type II and I models. Subsequently the Bianchi 
type IX case is discussed and the notion of an extended Kasner map for the Mixmaster singularity is introduced. 
Section ^ is devoted to discussing how the extended equations of this paper can be used to analyze the anisotropics of 
the CMB temperature. We end with a discussion in section |^ and mention further possible applications. In Appendix 
A we outline how the individual geodesies can be found if needed. 

In the paper, latin indices, a, b, c, ... — 0, 1, 2, 4 denote spacetime indices while greek indices, a, {3, ... — 1, 2, 3 denote 
spatial indices in the orthonormal frame. 

2. EXTENDED ORTHONORMAL FRAME APPROACH 

In this section we derive the extended system of first-order differential equations that governs the evolution of 
SH universes and their geodesies. We introduce a group- invariant frame {go} such that eg = n is the unit normal 
vector field of the SH hypersurfaces. The spatial frame vector fields are then tangent to these hypersurfaces. The 
gravitational variables are the commutation functions of the orthonormal frame, which are customarily labeled 

{H,(7af},^a,napTaa} , (1) 

(see WE, equation (1.63)). The Hubble scalar H describes the overall expansion of the model, (Ta/3 is the shear 
tensor and describes the anisotropy of the expansion, naf3,aa describe the curvature of the SH hypersurfaces, and Qa 
describes the angular velocity of the frame. The evolution equations for these variables are given in WE (equations 
(1.90)-(1.98)). To be able to incorporate a variety of sources, we use the standard decomposition of the energy- 
momentum tensor Tab with respect to the vector field n. 

Tab = tmaUb + 2qi^anb) + P {gab + naUb) + TTab , (2) 

where 

qaU" = , TTab-n^ = TTa" = , 1T[ab] = . (3) 

Hence, relative to the group invariant frame, we also have the following source variables 

{t^,P,qa,T^ati} ■ (4) 

We now normaliz^ the gravitational field variables and the matter variables with the Hubble scalar H . We write: 

}/H, (5) 



^See WE, page 112, for the motivation for this normalization. 
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and 

{Q.,P,Qa,Tlaf3} ^ {^J',P,qa,^^af3} /iH^ ■ (6) 

These new variables are dimensionless and are referred to as expansion- normalized variables. By introducing a new 
dimensionless time variable r according to 

the equation for H decouples, and can be written as 

H' = -{l + q)H, (? = 2S2 + i(r! + 3P) , (8) 
where a prime denotes differentiation with respect to r. The scalar E is the dimensionless shear scalar, defined by 

and q is the deceleration parameter of the normal congruence of the SH hypersurfaces^. The evolution equations for 
the dimensionless gravitational field variables follow from equations (1.90)-(1.98) in WE, using (^-(^). 
We will now consider the geodesic equations, 

fc-^Vafc' = , (10) 

where k'^ is the tangent vector field of the geodesies^ We can regard an individual geodesic as a curve in a spatially 
homogeneous congruence of geodesies, in which case the orthonormal frame components of its tangent vector field 
satisfy 

ea(r) = 0. (11) 

We now use equations (1.15) and (1.59)-(1.62) in WE to write ( [lO| ) and ( pT| ) in the orthonormal frame formalism, 
obtaining 

k"k° = ^aa^pk^k^ - H {k"f , (12a) 
fcOfc" = -fc/^ (cr"^ + HS°'p) fc° + e°'fsun^k^k^' - apkl^k" + a" (fc^fc") , (12b) 

where an overdot denotes differentiation with respect to t, the cosmological clock time (synchronous time). We now 
introduce energy-normalized geodesic variables 

k" 

A- = - , (13) 

where £ = k'^ is the particle energy. The vector K" satisfies KaK" = 1 for null geodesies, < 1 for timelike geodesies, 
and > 1 for spacelike geodesies. For null geodesi es, the variables Ka correspond to the direction cosines of the 
geodesic. The equation for the energy £, equation (12a) decouples and can be written as 

£" = -(l + s)£, (14) 

where 

S ^ -1 + KaK" + T.afiK"K'^ . (15) 

We now summarize the extended system of equations in dimensionless form. 



*The equation for q generalizes equation (5.20) in WE. 

^For many purposes in SH cosmology, it is sufficient to consider only the geodesic tangent vectors, and not the coordinate 
representation of the geodesies themselves. If specific coordinates are introduced, the geodesies can be found by the methods 
outlined in appendix A. 
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Evolution equations 



Kf3 - -(2 - 9)S„/3 + 2e''%j:p)^,R, + , (16a) 

= qN^p + 2S(/7V^)^ + 2e''%Npj^R, , (16b) 

= qA„ - E/A^ + e/^A^i?, , (16c) 

i^;, = {s + ApRP) - ^^pK" - to^p^R^K^ - e^t3^N\K'K^ - {KpK^) , (16d) 

Constraint equations 

n = l-T? -K , (16e) 

3Q„ = 3S/A^ - e/^E/TV^, , (16f) 

= iVa'^A/J , (16g) 

where the spatial curvature is given by 

'5„/3 = B„;3 - i (S/) - 2e^%7V;3)^^, , (17) 

= , (18) 

with 

= 27V/7V,,;5 - (TV/) 7V„^ . (19) 

Accompanying the above system of equations are, if necessary, equations for matter variables. For example, if the 
source were a tilted perfect fluid, additional equations for the tilted fluid 4-velocity would have to be added. 

Note that the null geodesies, characterized by KaK" = 1, define an invariant subset. This is easily seen from the 
auxiliary equation for the length of the vector K", 

{K^K°')' ^~2{1- K^K°'){l + s) . (20) 

From now on we will restrict our considerations to null geodesies, in which case the expression for s simplifies to 

s = T^apK'^K'^ . (21) 

A. Examples: Some non-tilted perfect-fluid models 

For non-tilted perfect fluid models, the 4-velocity of the fluid, u, coincides with the normal vector field n and 
Qa = Ilab = 0. It will also be assumed that the cosmological fluid satisfies a linear barotropic scale-invariant equation 
of state, p = (7 — l)/i, or equivalently, P = (7 — where 7 is a constant. From a physical point of view, the most 
important values are 7 = 1 (dust) and 7 = | (radiation). The value 7 = corresponds to a cosmological constant 
and the value 7 = 2 to a "stiff fluid". Here it is assumed that < 7 < 2. Our focus will be on diagonal Bianchi 
models. These are the class A models, and the Na" — models of class B, i.e. models of type V and special models 
of type Ylh (see Elhs & MacCallum @). 



Class A models 



For the class A models, Aa — 0, it is possible to choose a frame such that Nap = diag(A^i, -/V3), Ra ~ 0, and 

Ea/3 = diag(S+ + V3S_ , E+ - VSE- , -2E+) , (22) 

(see WE, page 123). Here we have chosen to adapt the decomposition of the trace-free shear tensor Ec./? to the third 
direction, rather than the first direction, as in WE. The anisotropic spatial curvature tensor ^Sap is also diagonal and 
we label its components in an analogous way: 

^Sap = diag (S+ + , S+ - , -25+) . (23) 

With the above choice of frame, (Iq) leads to an extended system of equations of the form: 
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Evolution equations 



where 



S]'±--(2-g)S]±-5± , (24a) 

N[^{q + 2E+ + 2V3Y.-)Ni , (24b) 

N^^{q + 2E+ - 2 V3I]_ )iV2 , (24c) 

A^^ = (q-4E+)7V3 , (24d) 

K[ = {s - S+ - V3E_)Xi + (A^2 - N3)K2K3 , (24e) 

i^^ = (s - S+ + \/3E-)X2 + (A^3 - Ni)KiKz , (24f) 

ir^ = (s + 2S+)i^3 + (A^i-A^2)i^ii^2 , (24g) 



q = i(37 - 2)(1 ^K) + 1(2 - 7)(S^ + , (25) 

s = (1 - 3i^|)E+ + V3 (i^i^ - if2^) S_ , (26) 

5+ = i [(7Vi-iV2)2-7V3(27V3-iVi-7V2)] , (27) 

S- = ^{N2-Ni){N^-Ni-N2) , (28) 

X = [Aff + Nl + iV| - 2 (Ari7V2 + 7V2iV3 + TVaA^i)] • (29) 
The density parameter is defined by 

17 = 1 - - si - X . (30) 



Diagonal class B models 



For the non-exceptional class B models with rioP = (denoted Ba and Bbi in Ellis & MacCallum pages 115,121- 
122), we can choose the spatial frame vectors Bq, so that the shear tensor Hap is diagonal, Ra — — (0,0, As), 

and the only non-zero compo nent s of Ngg a re A'^12 = A^2i- These models correspond to Bianchi type V and special 
type Wh models. Equations (16b) and (16c) imply that (-/Vi2/^3)' = 0, i.e. we can write 



where h is the usual class B group parameter. For convenience, we introduce a new parameter k according to 

1 



k = 



(31) 



(32) 



VI -3/i 

The type V models are ch arac terized by /c = 0, while fc = 1 corresponds to type VIq models, which are actually of 
Bianchi class A. Equation (16f) leads to restrictions on the shear tensor Eq/j, which can be written as 



= diag (-fc + V3VT~fc2 , -fc - \/3VT~~fc2 , 2fc)s> 
We now introduce a new variable A, and rewrite (pl|) in terms of k, obtaining 



Using 



A3 = Vl - k^A , N12 = VSkA . 
and (]3^), the extended system ( p^ reduces to the following set: 

Evolution equations 



SI 



-(2- g)Sx - 2/1:^2 



A' = (g + 2fcEx)^ 



K'2 



[k - V3VT~fc2)Ex + ( v/l-P + Vsfc) AK; 



(fc + VsVl - fc^)Sx + ( VT~~fc^- VSfc) AK: 



K2 



^/^ - k'^AK3 + 2kT., 



V3{Kf -~Kl)k+ ^1-k^ 



A 



(33) 



(34) 



(35a) 
(35b) 

(35c) 
(35d) 
(35e) 
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where 

g = 1(2 -7)SU 5(37 -2)(1-A2), (36) 



s= (fc - V3 - k^)Kf + {k + V3 - k^)Ki + 2kKl Ex- (37) 
The density parameter is given by 

n = l-Y.l- . (38) 

3. STRUCTURE OF THE EXTENDED SYSTEM OF EQUATIONS 

We now give an overview of the structure of the combined system of gravitational and geodesic equations. For 
simphcity, we only consider the non-tilted perfect fluid models described in section ^j. The basic dimensionless 
variables are 

X = {^^p , N^p , A^} , (39) 
K = {K^} . (40) 

We have shown that the Einstein field equations lead to an autonomous system of differential equations of the form 

x' = f(x), (41) 

(see (p4^-d) and (|3^a-b)). The geodesic equations lead to an autonomous system of differential equations of the form 

K' = h(x,K) , (42) 

which is coupled to ( pl| ) (see (p^-g) and (|35|c-e)). The geodesic variables Ka also satisfy the constraint 

K^K" = 1 , (43) 

and hence define a 2-sphere, which we will call the null sphere. In the context of cosmological observations, one can 
identify the null sphere with the celestial 2-sphere. We will refer to the entire set, equations (^)-(^) for x and K, 
as the extended scale- invariant system of evolution equations, or briefly, the extended system of equations. 

There are also two variables with dimension, namely the Hubble scalar H and the particle energy £. These scalars 



satisfy the decoupled equations (|^) and (14). They are thus determined by quadrature once a solution of the extended 
system of equations has been found. 

We now discuss the structure of the state space of the extended system of equations (^T|)-(^2|). The fact that the 
gravitational field equations ( ^l| ) are independent of K implies that the state space has a product structure, as follows. 
For models of a particular Bianchi type the gravitational variables x belong to a subset B of R" (see WE, section 
6.1.2 for Bianchi models of class A). Because of the constraint (^), the extended state space is the Cartesian product 
B X S^, where S'^ is the 2-sphere. The orbits in B lead to a decomposition of the extended state space into a family of 
invariant sets of the form {F} x S'^, where F is an orbit in B. Given a cosmological model U, its evolution is described 
by an orbit Fj/ in B. The orbits in the invariant set {Tu} x S'^ then describe the evolution of the model and all of its 
null geodesies. We shall refer to {Tu} x S"^ as the geodesic submanifold of the model U in the extended state space 
B X S^. In physical terms, with the null sphere representing the celestial sky, the geodesic submanifold of a model U 
determines the anisotropy pattern of the CMB temperature in the model U (see section ||). 

An advantage of using a scale-invariant formulation of the gravitational evolution equations is that models admitting 
an additional homothetic vector field, the so-called self-similar models, appear as equilibrium points (see WE, page 
119). The equilibrium points of the field equations are constant vectors x = Xq satisfying f (xq) = 0, where f is the 
function in (|4l|). In this case, the geodesic equations, 

K' = h(xo,K), (44) 

form an independent autonomous system of differential equations. The equilibrium points of the extended system 
(^-(^^ are points (xo,Ko) that satisfy 

f(xo) = 0, h(xo,Ko) = 0. (45) 

Knowing the equilibrium points of the field equations (see WE, section 6.2 for the class A models) one simply has to 
find the equilibrium points of the geodesic equations in (Q). The fixed point theorem for the sphere guarantees that 
the system of geodesic equations for self-similar models has at least one equilibrium point on the null sphere. Since 
the null sphere can be identified with the celestial 2-sphere, equilibrium points of the extended system of equations 
correspond to the existence of geodesies in fixed directions, i.e. purely "radial" geodesies. 
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4. EXAMPLES OF EXTENDED DYNAMICS 



In this section we will consider some examples of self-similar and non-self-similar models. For self-similar models, 
the extended system of equations reduces to (|^, and it is possible to visualize the dynamics of the geodesies. The 
most important self-similar models are those of Bianchi type I and II, namely the flat Friedmann-Lemaitre model, the 
Kasner models and the Collins-Stewart LRS type II model (see Collins & Stewart ||]), since these models influence 
the evolution of models of more general Bianchi types. For non-self-similar models, the dimension of the extended 
system of equations is usually too large to permit a complete visualization of the dynamics although one can apply 
the standard techniques from the theory of dynamical systems. In the simplest SH cases, however, one can visualize 
the dynamics, and as an example of non-self-similar extended dynamics, we will consider the Bianchi type II LRS 
models. We will end the section with a discussion of the Bianchi type IX models. 



A. Self-similar models 

The flat Friedmann-Lemaitre model 

The flat FL model corresponds to the following invariant subset of the extended system of equations for class A 
models: S+ — S_ = A^i = A^2 = -^3 = 0. The remaining equations in ( p^ ) are just 

K^O, a ^1,2, 3. (46) 

Thus, all orbits corresponding to null geodesies are equilibrium points and the null sphere is an equilibrium set. This 
fact implies that all null geodesies are radial geodesies. 



Kasner models 



Although these are vacuum models, they are extremely important since they are asymptotic states for many of 
the more general non- vacuum models. The models correspond to the Bianchi type I invariant vacuum subset of the 
extended system of equations for the class A models: A^i = A^2 = -^3 = ,1]"^ + = 1, where and I]_ are 
constants. 



The remaining equations are the geodesic equations (24e) - (241) with iVi = A^2 = A3 = and with s given by 



(Pq). We note that these equations are invariant under the discrete transformations 

iKi,K2,K3) ^ {±Ki,±K2,±K3) . (47) 

The constant values of T.^ and T._ determine the so-called Kasner parameters according to (see WE, equation 
(6.16) with 1,2,3 relabeled as 3,1,2) 

= i f 1 + S+ ± , P3 - i (1 - 2S+) . (48) 



Pi, 2 = 

One can also label the Kasner solutions using an angle (p, defined by E+ = cos(^ and I]_ = sin(^. All distinct models 
are obtained when assumes the values < (/s < ^. The equilibrium points for these equations are listed in table 1, 
together with their eigenvalues. In the LRS cases pi = P2 7^ there is a circle C12 of equilibrium points. 

For each set of Kasner parameters pi,P2, and ps the geodesic equations admit local sinks and local sources, which 
can be identified by considering the signs of the eigenvalues in table 1. It turns out that these local sinks/sources 
are in fact global, i.e. attract/repel all orbits, and hence define the future/past attractor. The reason for this is the 
existence of monotone functions that force all orbits to approach the local sinks/sources into the future/past. For 
example, for models with S+ 7^ we have the function 

A1A2 

The future and past attractors are listed in table 2 for the three cases pi ^ p2 ^ {i.e. ip — 0), pi > p2 > > P3 
(i.e. < If < ^) and P2 = P3 = {i.e. </? = f). 

In figure |l] we show the orbits corresponding to null geodesies in the Kasner models for the three cases Pi = P2 ^ 
{ip = 0), pi > p2 > > p3, and P2 = P3 = ((/s = ^). Due to symmetry, it is sufficient to show the subset of the null 
sphere defined by i^i, ^^2, ^^3 > 0. 
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Eq. point 




Eigenvalues 




A" 


(±1,0,0) 


-3(p2 -Pi) , -3(P3 -Pi) 


Pl^P2^ P3 


— 

p 

^2 


(0,±1,0) 


-3(pa -P2) , -3(pi -p2) 




P.,+ 
J 3± 


/'n n ^1 ^ 

(U, U, ±ij 


— 3(p2 — pa) , — 3(pi — pa) 


Pi = P2 / 


Cl2 


(cos ^, sin 0) 


, 3 


^a 


(0,0, ±1) 


-3 , -3 


Pi = P2 = 


Cl2 


(cos ^, sin i/), 0) 


, -3 


^a 


(0,0, ±1) 


3 , 3 


Table 1. The equilibrium points of the geodesic equations for Kasner models, written in terms of pi, p2, and ps as defined 
in the text. The parameter is a constant satisfying < i' < 2n. The eigenvalues for other LRS models than pi — p2 can be 
found by appropriate permutations. 


Kasner parameters 




Past attractor 


Future attractor 


Pi > P2 > > P3 






^Pa^UP3-} 


Pi = P2 / 




Cl2 




P2 = Pa = 






C23 



Table 2. The past and future attractors in the state space for null geodesies in the Kasner models whose parameters satisfy 
pi > P2 > > Pa. The results for other ordering of parameters can be obtained by appropriate permutations. 



The Collins-Stewart LRS type II solution 

The Collins-Stewart model corresponds to the following subnianifold^ of the extended system of equations: 

E+ - i(37 - 2) , S_ = , = iV2 = , iVg = | V(2 - 7)(37 " 2) , (50) 

with I < 7 < 2. Due to the symmetries, we need only consider K3 > 0. The equilibriimr points and sets are listed in 
table 3. 

The equilibrium set C12 is the source, while the the equilibrium point P3 is a stable focus. Note that is an 
increasing monotone function. The dynamics of the null geodesies is shown in figure ^ Note that there are no changes 
in the stability of the equilibrium points for | < 7 < 2. 

B. Non— self— similar models 

The previous examples are simple in the sense that we only had to consider the geodesic part of the extended system 
of equations. For non-self-similar models, the full system has to be considered, which means that the dynamics will 
in general be difficult to visualize due to the high dimensions of the extended state space. To illustrate the ideas, we 
consider the null geodesies in Bianchi type I and II LRS models. The behavior of geodesies in the Mixmaster model 
is also discussed. 



Note the incorrect numerical factor on page 131 in WE. 



Eq. point 


Ki,K2,K3 


Eigenvalues 


Pa 


0,0,1 


f(2-7),-^(37-2±2ib) 


C12 


costp, sintp, 


f(2 + 7),f(37-2),0 



Table 3. The equilibrium points and sets for the null geodesic equations in the Collins- Stewart LRS type II solution. The 
parameter 7/) satisfies < < 2n. Note that two of the eigenvalues for the equilibrium point C12 are complex. The constant b 
is given by b = ^(87 - 2)(2 - 7). 
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Figure 2. The dynamics of the null sphere for the Collins-Stewart LRS type II model, as viewed from the positive 3-axis. 



LRS Bianchi type I and 11 models 

The type II LRS models correspond to the invariant subset S_ = 0, iVi = A^2 = of the extended system of 
equations (^J) for the class A models, while the type I models, in addition, require A^3 = 0. For null geodesies, the 
extended system is five dimensional (four for type I), with one constraint KaK°' — 1. Defining 



Ki — Rcosx , K2 — Rsiii-x , (51) 



where 

leads to a decoupling of the x-equation, 
leaving a reduced extended system 



with 



R = , (52) 

X' = N^K^i , (53) 

SV--(2-g)I]+ + i7V| , (54a) 
= (g - 4E+)iV3 , (54b) 
K's = SJ^+il - Kl)K3 , (54c) 

<? = i(37 - 2)(1 - ^Ni) + 1(2 - . (55) 

The state space associated with (^) is the product set B x [0, 1], where B is the state space of the Bianchi type II 
LRS cosmologies (or type I, in the case — 0), associated with the subsystem of (|5^a)-(|54[3). In this representation 
the null sphere is replaced by the single geodesic variable K3, with < K3 < 1. The remaing two geodesic variables 
are given by (|5^). We refrain from giving the various equilibrium points and their eigenvalues. Instead we give the 
three-dimensional extended state space of ( p4[ ) in figure and the two-dimensional invariant set = in figure 
In figure |^ we have simply shown the skeleton of the state space, i.e. the equilibrium points and the various 
heteroclinic orbits that join the equilibrium points. The figures depict the situation when | < 7 < 2 since there are 
no bifurcations for this interval. The sources and sinks can be deduced from the figures. A detailed picture of the 
orbits in the gravitational state space K3 = is given in WE (see figure 6.5). We note that the orbits in the invariant 
set K3 = 1 are identical to those with K3 = 0. Knowing the orbits in — and — 1, one can visualize the 
structure of the geodesic submanifolds - they are vertical surfaces of the form {F} x [0, 1], where F is an orbit in the 
subset K3 = 0. 
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I N3 



(a) (b) 

Figure 3. The dynamics of the extended system of equations for (a) Bianchi type I LRS models and (b) Bianchi type II LRS 
models. 

Comments on Bianchi type IX models 

It was recognized a long time ago that the oscillatory approach to the past or future singularity of Bianchi IX 
vacuum models, the so-called Mixmaster attractor, displays random features, see e.g. Belinskii et al. and hence is 
a potential source of chaos. This behavior is also expected in non- vacuum Bianchi models with various matter sources 
(see section 6.4.1 in WE and references therein). Numerical studies of the governing equations of vacuum Bianchi IX 
models toward the initial singularity have shown that the variables S]± and the Na remain bounded. These studies 
have also shown that the projection of the orbits onto the E-t-plane is given, at least to a high accuracy, by the Kasner 
map (see WE, section 11.4.2). The transition between two different Kasner states is described by a vacuum Bianchi 
type II orbit except when the Kasner state is close to an LRS Kasner model where this approximation is no longer 
valid. When discussing the Mixmaster attractor, one is usually discussing individual orbits. Thus a corresponding 
discussion for the extended system implies a discussion about the Bianchi type IX geodesic submanifold. Precisely as 
an individual Bianchi type IX orbit can be approximated by a sequence of Bianchi type II orbits, one can approximate 
a type IX geodesic submanifold with a sequence of type II geodesic submanifolds. The stable equilibrium points within 
the type II geodesic submanifold reside in the type I geodesic boundary submanifold of these models and correspond 
to geodesies in the 1,2 or 3 directions, modulo sign, depending on the particular Kasner point. 

We will only consider such sequences of Kasner states for which the Kasner models are not close to any LRS models. 
As the evolution progresses, the r-time that the system spends close to a Kasner state, a so-called Kasner epoch, 
becomes successively longer and should thus be well described by the appropriate equilibrium point. If we assume 
that during a certain Kasner epoch the qualitative behavior of a geodesic is given by the stable equilibrium point of 
the extended system of equations for these models, we can extend the Kasner map to include the stable direction of 
the geodesic. Since we are excluding the LRS Kasner models there will never appear any equilibrium sets as they 
only when cp is a multiple of 7r/3. The direction of stability, modulo sign, is given in table 4 as a function of ip. 
These results follow from the general stability of the equilibrium points given in table 1, by changing the signs of the 
eigenvalues since the models are approaching the the initial singularity, i.e. t — s- — 00. 

Starting with a geodesic whose tangent vector satisfies Ki,K2,K^ > in a given Kasner epoch with < (pini < tt/S, 
the stable geodesic direction is the 1-direction. The system then evolves, according to the Kasner map, into a state 
with ^ < ipfin < TT. Depending on the initial value (p-mi, the stable geodesic direction can either stay the same 
(arccos (-jl) < ipmi < f ) or change to the 2-direction (0 < (^i„i < arccos (3^)). This process is then repeated as the 
state changes again. This extended Kasner map is shown in figure ^. In the figure, a whole sequence of Kasner states 
is also shown where the stable geodesic directions are given by the sequence 1, 1, 2, 2, 1, 3. The above discussion of the 
behavior of geodesies toward the Mixmaster singularity is based on the assumption that as the system changes from 
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Range of Lp 


Stable geodesic direction 


< < 27r/3 


1 


27r/3 < 1/3 < 47r/3 


3 


47r/3 < (/3 < 27r 


2 



Table 4. The stable geodesic direction for the geodesies for a particular Kasner epoch for different (^'s. By assumption we 
exclude all the LRS Kasner models, i.e. models when ip is, a. multiple of 7r/3. 




Figure 4. The extended Kasner map including the stable geodesic directions. The sequence of Kasner states depicted has 
the following sequence of stable geodesic directions: 1,1,2,2,1,3. 



one Kasner epoch to another, the geodesies are not affected. This means that a tangent vector to the geodesic with 
^1,2,3 > can never evolve into a tangent vector with one or more of the K^s negative. 

This assumption is rather crude since the change of Kasner epochs is approximately described by a vacuum Bianchi 
type II orbit, for which the geodesic equations, if viewed as separate from the field equations, are non-autonomous. 
Taking this into account limits the predictability of the "extended Kasner map" in that it cannot predict if a geodesic 
evolves in the positive or negative direction of the stable geodesic direction. We also note that it is only in r-time that 
the system spends longer and longer time in each Kasner epoch. In synchronous time, the interval becomes shorter 
and shorter. 

From the above discussion, it is expected that there will be some kind of geodesic chaos in the development toward 
the initial and final singularity. To substantiate this, one would need, in addition to further analytical results, careful 
numerical studies. We believe that the extended system of equations, as presented in this paper, may be very well 
suited for such an analysis. The next step would be to study the extended system of equations for Bianchi type II 
vacuum models. 



5. TEMPERATURE DISTRIBUTION 

In this section we describe how the extended system of equations (|4l|)-(|4^), together with the decoupled energy 
equation (|T^), can be used to study the temperature of the CMB in an SH universe. We regard the photons of the 
CMB as a test fluid, i.e. one which is not a source of the gravitational field. It is possible to include the effect of 
the CMB photons on the gravitational field by considering two non-interacting fluids, radiation and dust, using the 
approach of Coley & Wainwright |10|. We will not do this since the effects of the radiation fluid is not expected to 
change our results significantly. To obtain the present temperature of the CMB, the photon energies are integrated 
along the null geodesies connecting points of emission on the surface of last scattering to the event of observation at 
the present time. To simplify the discussion, it is assumed that the decoupling of matter and radiation takes place 
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instantaneously at the surface of last scattering. The matter of the background cosmological model is assumed to be 
described by dust, i.e. p = 0. 

By the following simple argument we can approximate the interval of dimensionless time, At, that has elapsed from 
the event of last scattering until now. If the radiation is thermally distributed, its energy density /ir, as derived from 
the quantum statistical mechanics of massless particles, satisfies /i^ cx where T is the temperature of the radiation 
(see Wald [0, page 108). A non-tilted radiation fluid satisfies fi^- cx exp(— 4t), which implies 



To 



-At 



(56) 



Here To and Te are the temperature at the present time and at the surface of last scattering respectively. Assuming 
that the process of last scattering took place when Tg « 3000 K, and that the mean temperature of the CMB today 
is To ~ 3, it follows that At ^ 7. This corresponds to a redshift of about z « 1100. 

The temperature of the CMB can now be found as follows. Introduce a future-pointing null vector k which is 
tangent to a light ray at a point on the CMB sky. The current observed temperature To of the CMB is given by (see, 
for example, Collins & Hawking h], page 313 ) 



To 



From Eqs. (|T|) and (|l]) it follows that 

ToSxp 



To 



jugk'')^ ^ £(to)_ 
(uak-)^ £{To) 



^,,p{T)K^{T)KP{T))dT 



(57) 



(58) 



This formula gives the temperature at time t = tq in the direction specified by the direction cosines Ka{TQ). We 
introduce angles 6, ip by 



Ki (to ) = sin 9 cos ip 
K2 {Tf)) — sin 9 sin p . 



(59) 



to describe positions on the celestial sphere. Note that to obtain a correspondence with the spherical angles defining 
the direction in which an observer measures the temperature of the CMB, one has to make the transformation 
9—i-TT — 9,ip—fip + TT. In this way, To is a function of the angles 9 and ip, i.e. 



To = T{9,p) , 



(60) 



which we call the temperature function of the CMB. 

The anisotropy in the CMB temperature can be described using multipole moments (see for example Bajtlik et al. 
H). The fluctuation of the CMB temperature over the celestial sphere is written as a spherical harmonic expansion. 



= ^ ^ aimYlm, {0, if) , 



1=1 m=-l 

where Tav is the mean temperature of the CMB sky. The coefficients aim a-re defined by 

AT 



T 



(61) 



(62) 



where * denotes complex conjugation, and the integral is taken over the 2-sphere (see for example Zwillinger [ p^ , 
pages 492-493). The multipole moments, describing the anisotropics in a coordinate independent way, are defined as 



1/2 



ai 



E 

\m=-l 



(63) 



The dipole, ai, is interpreted as describing the motion of the solar system with respect to the rest frame of the CMB. 
Therefore, the lowest multipole moment that describes true anisotropics of the CMB temperature is the quadrupole 
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moment, 02. Current observations provide an estimate for 02 as well as for the octupole moment 03 (see Stoeger et 
al. H). 

In order to compute T{9, ip) and the multipole moments 02 and 03 for a particular cosmological model, one has to 
specify the dimensionless state, x(to) — Xq, of the model at the time of observation, Tq, and the direction of reception 
Ka{To), which determines the angles 9 and (p on the celestial sphere via (|59|). The solution x(r), Kair) of the 
extended system of equations (^ ) and (^2|), determined by the initial conditions x(ro) and Ka{To), is substituted in 
( ^ ) , which determines the temperature function T {9, (p) . The multipoles 02 and 03 are then calculated by integrating 
over the 2-sphere (see ( |6^ ) and (|6^)). In this way the multipole moments can be viewed as functions defined on the 
dimensionless gravitational state space, with the time elapsed At since last scattering as an additional parameter: 

0,2 = 02 (xo; At) , 03 = 03 (xo; At) . (64) 

The extended equations can be used in three ways to obtain information about T{9, if) and the multipoles 02 and 03, 
as follows. 

i) Apply dynamical systems methods to the extended equations to obtain qualitative information about the null 
geodesies and the shear, and hence about the temperature pattern of the CMB. 

ii) Linearize the extended equations about an FL model, and if possible solve them to obtain approximate analytical 
expressions for T(9^ ip), a-i and a3, which are then valid for E << 1. 

iii) Use the full non-linear extended equations to do numerical simulations, calculating T(0, (p), 02 and 03 for a given 
point Xo in the gravitational state space, not necessarily satisfying So << 1. As described above, one calculates 
the value of T (0, tp) at each point of a grid covering the celestial sphere and then integrates numerically over 
the sphere to obtain 02 and a3. 

One can use the observational bounds on 02 and 03, together with the results of ii) and iii) above, to determine bounds 
on the anisotropy parameters associated with the shear and the Weyl curvature, So and Wo- A necessary condition 
for the model to be close to FL at the time of observation is that So and Wo are small (see Nilsson et al. |lj]). 

We are in the process of applying the above methods to SH models of various Bianchi group types. Preliminary 
results on Bianchi VIIq models are given in Nilsson et al. |Q , where it is shown that the observational bounds on 02 
and 03 do not necessarily imply that the Weyl parameter Wo is small. Thus, an almost isotropic CMB temperature 
does not imply an almost isotropic universe. An advantage of the above methods is that they are not restricted to 
those Bianchi types that are admitted by the FL models. For example, we are studying the diagonal class B models of 
Bianchi type Vl/i as described by (p5[). The current state Xo can be described by So, k and Qq, and so the expression 
for the quadrupole has the general form 

a2 = a2(So, fc, rio; At) . (65) 

When one uses the method ii) and linearizes about the open FL model (So = k — 0), one obtains a formula of the 
form 

02 - ^ohi^o-, At) + kJ^in^, At) , (66) 

where I2 and J2 are expressed as integrals. This formula, which is valid for So << 1 and A; << 1, leads to bounds 
on the shear parameter So that are much weaker than those obtained in Bianchi types I and V (see Lim et al. [p5[). 
This result shows that the bounds obtained for the anisotropy parameters in Bianchi type 1 and V models (see, for 
example Bajtlik et al. Q), which seem to have been taken for granted as being typical, are misleading. 

6. DISCUSSION 

In this paper we have shown that in the case of spatially homogeneous models, the field equations can be augmented 
with the geodesic equations, producing an extended set of first-order evolution equations whose solutions describe not 
only the evolving geometry but also the structure of the geodesies. Examples of the dynamics of geodesies in some 
self-similar models, and in a simple non-self-similar model, in order to show the predictive power of the approach, 
were given. We also made some conjectures about the qualitative behavior of geodesies towards the initial Mixmaster 
singularity in Bianchi IX models by considering the geodesies structure of the Kasner models. In light of the chaotic 
nature of the Mixmaster singularity, it is expected that the geodesies will also have some sort of chaotic behavior. 
To confirm these speculations, we point out that further analytic studies and thorough numerical investigations are 



13 



needed. We believe that the formulation of the combined field equations and geodesic equations of this paper is very 
well suited for this purpose. 

The most physically interesting aspect of the extended set of equations is the possibility of shedding light on one 
of the fundamental questions concerning the CMB, namely, does a highly isotropic CMB temperature imply that the 
universe can accurately be described by an FL model? As mentioned in section ^, the preliminary indications are 
that the situation is less clear-cut than was previously thought. 

We now conclude by listing some related research topics. One can easily generalize the present formalism to spatially 
self-similar models and it is also easy to include other sources, e.g. , two non-interacting fluids, a cosmological constant, 
or magnetic fields. It should also be possible to extend the system to include polarization. It would also be interesting 
if the current formulation could be generalized to facilitate a study of the temperature of the CMB for density 
perturbed models or models with closed topologies. 



APPENDIX A: OBTAINING THE GEODESICS 



For many physical purposes in SH cosmology, it is sufficient to know the tangent vector field k of the geodesies. 
Nevertheless, it might be interesting to find the geodesies themselves. To do so, spacetime coordinates x'^ (/i = 0, 1, 2, 3) 
must be introduced. Once this is done, a geodesic can locally be described as a curve a;^(A), where A is the affine 
parameter of the geodesic. With an SH geometry, it is natural to adapt the coordinates to the structure imposed by 
the SH condition. Expressed in terms of coordinates, the orthonormal frame can then be written as 

eo = iV(0-i(ai-iV'(<,x)E,(x)) , = e„'(t)E,(x) , (Al) 

where N and N = iVE, are the lapse function and the shift vector field respectively (for restrictions on the shift vector 
field, see Jantzen & Uggla @). The spatial frame {Ei(x)}, where i = 1,2,3 (in this appendix, and only here, latin 
indices i,j,k... = 1,2,3 denote spatially homogeneous time independent frame indices since it is only here that this 
type of frame is used), tangent to each hypersurface, is not only invariant under the action of the Bianchi symmetry 
group but has structure or commutator functions C\j which are constants throughout the spacetime, defined by 

[E„E,] =C% Efc . (A2) 

It is possible to construct the orthonormal frame explicitly in terms of local coordinates {t, x'} {i = 1, 2, 3) adapted to 
the SH hypersurfaces. The spatial frame {E^} is characterized by the Lie dragging condition Ce^Ei = which implies 
the time independent local coordinate expression E^ = Ei^ {x^) dj for the invariant spatial frame. Explicit coordinate 
expressions for Ei^ {x^ ) follow from the representation of the left invariant vector fields in canonical coordinates of the 
second kind in Jantzen [|l^,^. The relation between the orthonormal frame components of the tangent vector field 
to the geodesies and its corresponding coordinate components yields the relations 

{H£)-^ , (A3a) 

^[-iV^(T,x) + ff-i(^X"(r)e„^(T))]i?/(x) , (A3b) 

j OL 

where the K°' are the energy-normalized components of the geodesic tangent vector field. Here one can choose an 
automorphism adapted shift in order to set certain components of {t) to zero, see Jantzen & Uggla ||l^. One 
usually sets the shift to zero. In this latter case one has to add equations governing some of the ea-' (r) components. 
These are obtained from the commutator relations and are consequences of the zero-shift gauge. We will not do this 
for the general case. Instead we will look at an example. 



dx' 
d7 



1. An example: Non-tilted class A models with zero shift 



In this case, with zero shift, equation (A3) becomes 

dX/dr = {HSy^ 



dx'/dr ^^H-'^K'^e~^"Ea'{x) 



(A4) 
(A5) 



14 



where eJ = 5^e The commutator equations (see Jantzen & Uggla I^6j or WE, chapter 10) yields 

df3i,2/dT = 1 + S+ ± V3I]_ , (A6) 
dPs/dr = 1 - 2E+ . (A7) 

The next step is to introduce H^^e~^" as new variables, but note that they are not independent of (see WE, 
chapter 10). This results in the complete extended system of equations, which is necessary in order to obtain the 
individual geodesies. 
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